//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   DD_2D_Equations.i.hh
 * \author Jeremy Roberts
 * \date   Sep 8, 2011
 * \brief  
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//


#ifndef DD_2D_EQUATIONS_I_HH_
#define DD_2D_EQUATIONS_I_HH_

namespace slabtran
{

// Solve in a cell.
template <class MF, class PCF>
inline void DD_2D_Equations::solve(const int                  i,
                                   const int                  j,
                                   const double               sigma,
                                   const double               q,
                                   const Omega                &omega,
                                   const double               wt,
                                   const Moment_to_Discrete_t &M,
                                   const Space_Vector         &in_psi,
                                   MF                         &phi,
                                   Space_Vector               &out_psi,
                                   PCF                        &psi,
                                   const int                  angle,
                                   bool                       store_psi)
{
  // This implementation is nearly identical to the presentation
  // given in L&M eqs. 4-29 -- 4-33.

  // The abs(omega) lets us use the same code for both directions.
  double d_factor[] = {2.0 * fabs(omega[MU])  * d_mesh->inv_width(i, I),
                       2.0 * fabs(omega[ETA]) * d_mesh->inv_width(j, J)};

  // Calculate cell-centered angular flux.
  d_psi = (q + d_factor[I] * in_psi[I]
             + d_factor[J] * in_psi[J]) /
          (sigma + d_factor[I] + d_factor[J]);

  // If requested, store the cell-centered angular flux.
  if (store_psi) psi(i, j, 0, angle) = d_psi;

  // Calculate new exiting flux.
  out_psi[I] = 2.0 * d_psi - in_psi[I];
  out_psi[J] = 2.0 * d_psi - in_psi[J];

  // Loop through flux moments and add up the contributions.
  int cell = d_mesh->convert(i, j);
  int moment = 0;
  for (int l = 0, N = phi.legendre_order(); l <= N; l++)
    for (int m = -l; m <= l; m += 2, moment++)
      phi(cell, 0, moment) += d_psi * wt * M(angle, moment);
  /// \todo [PERFORMANCE] Use BLAS call here for loop (need direct access to
  ///       continuous data)
}

} // end namespace slabtran

#endif /* DD_2D_EQUATIONS_I_HH_ */

//---------------------------------------------------------------------------//
//              end of DD_2D_Equations.i.hh
//---------------------------------------------------------------------------//
